/*
 * kpm_ldos_kernels.cu
 *
 *  Created on: 2011-2-18
 *      Author: zhang
 */
#include <arguments.h>
#include <crs.h>
#include <stdlib.h>
#include <stdio.h>
#include <curand_kernel.h>
#include <pthread.h>
#include <sys/time.h>
#include <ell.h>
#include <math.h>

#define BLOCK_SIZE 128
#define BLOCKS 64
#define THREADS (BLOCKS*BLOCK_SIZE)


namespace LDOS
{
	/*********************
	 * firstly perform generate vector R1 = H * R,
	 * then perform v-v multiplication R*R1
	 */

	template<class T>
	__global__ void cukpmCreateR0(unsigned int r, T* R0, unsigned int size, T* Mu)
	{
		  int tid_global= blockDim.x * blockIdx.x + threadIdx.x;
		  unsigned int total_threads = gridDim.x * blockDim.x;

		  for(int i=tid_global;i<size; i+=total_threads){
			  R0[i] = 0;
			  if(i == r){
				  R0[r] = 1;
				  Mu[0] = 1;
			  }
		  }
	}

	/**********
	 * H is in column storage.
	 */
	template<class T>
	__global__ void cukpmCreateR1(const ELL<T>* H, T* R0, T* R1, T* Mu)
	{
		  int tid_global= blockDim.x * blockIdx.x + threadIdx.x;
		  unsigned int total_threads = gridDim.x * blockDim.x;

		  for(int i=tid_global;i<H->rowsNum; i+=total_threads){
			T product = 0.0;
			for(int j=0;j<H->colsNum;j++){
				product += H->elements[j*H->rowsNum + i] * R0[H->cols[j*H->rowsNum + i]];
			}
			R1[i] = product;
		  }

	}

	template<class T>
	__global__ void cukpmCreateRn(const ELL<T>* H, T* R0, T* R1, T* Mu, unsigned int r, unsigned int n)
	{
		  int tid_global= blockDim.x * blockIdx.x + threadIdx.x;
		  unsigned int total_threads = gridDim.x * blockDim.x;

		  //perform the product of R1 * R
		  if(tid_global==0)
			  Mu[n-1] = R1[r];

		  for(int i=tid_global;i<H->rowsNum; i+=total_threads){
			T product = 0.0;
			for(int j=0;j<H->colsNum;j++){
				product += H->elements[j*H->rowsNum + i] * R1[H->cols[j*H->rowsNum + i]];
			}
			R0[i] = 2*product - R0[i];
		  }
	}


	template<class T>
	__inline__ __device__ T getKernelValue(int j, const int m){
		return (T)1/(m+1)*((m-j+1)*cos(PI*j/(m+1))+sin(PI*j/(m+1))*cos(PI/(m+1))/sin(PI/(m+1)));
	}

	template<class T>
	__inline__ __device__ T cukpmCreateDOSValue(const T* h_U, const int m, const T x){
		int j;
		T Re=0;
		T T0,T1;

		j = 0;
		T0 = 1;

		Re += h_U[0] * getKernelValue<T>(0, m);//getKernelValue<T>(0, m);

		j = 1;
		T1=x;
		Re += 2 * h_U[j] * T1 * getKernelValue<T>(1, m);

		for(j=2;j<m;j++)
		{
		  T temp = 2 * x * T1 - T0;
		  T0=T1;
		  T1=temp;

		  Re +=2 * h_U[j] * T1 * getKernelValue<T>(j, m);
		}
		Re *= 1/(3.14*sqrt(1-x*x));
		return Re;
	}

	template<typename T>
	__global__ void cukpmCreateDOS(const T* h_U, const int m, const T* sample, T* values, const int s){
		  int tid_global= blockDim.x * blockIdx.x + threadIdx.x;
		  int total_threads = gridDim.x * blockDim.x;

		  for(int i=tid_global; i<s; i+=total_threads)
			  values[i] = cukpmCreateDOSValue<T>(h_U, m, sample[i]);
	}

	template<typename T>
	__global__ void cukpmCreateMeanDOS(const T* h_U, const int m, const T* sample, T* mean, T* geo_mean, const int s){
		  int tid_global= blockDim.x * blockIdx.x + threadIdx.x;
		  int total_threads = gridDim.x * blockDim.x;

		  for(int i=tid_global; i<s; i+=total_threads)
		  {
			  T value = cukpmCreateDOSValue<T>(h_U, m, sample[i]);
			  mean[i] += value;
			  geo_mean[i] += log(abs(value)+0.00000000001);
		  }
	}

	template<class T>
	__global__ void cukpmGetMean(T* mean, T* geo_mean, const unsigned int len, const int num){
		  int tid_global= blockDim.x * blockIdx.x + threadIdx.x;
		  int total_threads = gridDim.x * blockDim.x;

		  for(unsigned i = tid_global; i<len; i+=total_threads)
		  {
			  mean[i] /= (T)num;
			  geo_mean[i] = exp(geo_mean[i]/(T)(num-1));
		  }

	}

	template<class T>
	__global__ void cukpmCreateLastMu(T* d_U, const int m, const T* R1, const int r){
		if(blockDim.x * blockIdx.x + threadIdx.x == 1)
			d_U[m-1] = R1[r];
	}
}












